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ABSTRACT 

In  anticipation  of  the  AQN  201  nuclear  reactor's  operation  at  a 
power  of  1,000  watts  for  short  periods  of  time  and  at  20  watts  continuous- 
ly, a  prediction  of  the  behaviour  during  a  rapid  power  rise  is  obtained 
from  digital  computer  solutions  of  the  reactor  kinetic  equations.   The 
experimental  and  analytical  work  which  provided  the  coefficients  for  the 
equations  is  outlined.   Solutions  of  the  equations  indicate  the  peak  power 
which  may  be  achieved,  the  time  during  which  1,000  watts  may  be  obtained, 
and  the  temperature  rise  that  results. 
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1.   Introduction. 

The  AGN  201  nuclear  reactor  Is  a  low  power  reactor  designed  for 

educational  uses.   Usually,  it  is  licensed  to  operate  at  100  milliwatts 

and  serves  to  demonstrate  many  of  the  basic  characteristics  of  reactors „ 

Also,  at  this  power,  the  peak  neutron  flux  of  about  4.5  x  10  neutrons/ 

2 
cm  -sec  is  sufficient  for  the  creation  of  many  radioactive  isotopes  in 

practical  quantities. 

For  other  isotopes,  however,  it  is  desirable  or  necessary  to  have  a 
much  higher  neutron  flux  and  it  has  been  proposed  that  the  AGN  201  be 
operated  for  short  periods  at  a  power  of  1,000  watts.   Several  of  these 
reactors  have  been  operated  at  somewhat  higher  powers,  but  none  have 
operated  as  high  as  1,000  watts.   The  limiting  factor  on  elevated  power 
operation  is  the  lack  of  a  cooling  system.   With  no  means  of  increasing 
the  rate  of  heat  transfer  from  the  core,  the  overall  core  temperature  in- 
creases.  Since  the  reactor  has  a  negative  temperature  coefficient  of 
reactivity,  a  rise  in  temperature  reduces  the  available  excess  reactivity^ 
and  the  rate  at  which  power  may  be  further  increased  is  reduced. 

The  maximum  steady- state  power  is  that  power  at  which  the  rate  of 
heat  generation  in  the  core  is  equal  to  the  rate  of  heat  transfer  from 
the  core  and  the  corresponding  overall  temperature  of  the  core  has  re- 
duced the  initial  excess  reactivity  to  zero.   For  the  AGN  201,  maximum 
steady-state  power  has  been  estimated  to  be  approximately  20  watts. 
High  powers,  such  as  1,000  watts,  would  be  achieved  during  the  transient 
period  following  a  relatively  large  increase  in  reactivity  if  the  re- 
actor core  is  initially  at  room  temperature.   Under  such  conditions,  power 


can  be  increased  rapidly  by  several  orders  of  magnitude  before  the 
temperature  has  risen  sufficiently  to  make  the  reactor  sub-critical. 

The  transient  behaviour  of  a  reactor  can  be  predicted  using  reactor 
kinetic  equations  based  on  various  approximations.   In  general,  the  equa- 
tions cannot  be  solved  analytically  due  to  the  coupling  between  them, 
and  analog  computer  solutions  are  frequently  obtained  for  cases  where 
variations  of  the  dependent  variables  are  not  large.   When  the  dependent 
variables  change  by  orders  of  magnitude,  numerical  methods  may  be  ap- 
plied to  solve  the  equations  using  a  digital  computer. 

It  is  the  purpose  of  this  investigation  to  establish  a  set  of  re- 
actor kinetic  equations  for  the  AGN  201  and  obtain  a  solution  to  the 
equations  using  a  digital  computer  for  the  case  of  a  rapid  rise  in 
power.   Features  of  significance  to  be  obtained  from  the  solution  in- 
clude an  estimate  of  the  peak  power  that  may  be  achieved  (or  in  the  run- 
away accident  sense,  a  peak  power  which  will  not  be  exceeded),  the  length 
of  time  that  a  power  of  1,000  watts  may  be  maintained,  and  the  tempera- 
tures which  would  result.  The  latter  is  significant  since  there  is  a 
thermal  fuse  in  the  core  which  allows  the  bottom  half  of  the  core  to  drop 
away  in  the  event  of  overheating.   The  fuse  is  always  at  a  higher  tempera- 
ture than  the  rest  of  the  core. 

The  coefficients  of  the  kinetic  equations  are  obtained  from  experiment- 
al data  or  from  calculations.   The  available  excess  reactivity  of  the  re- 
actor, the  rate  at  which  reactivity  can  be  added,  and  an  estimate  of  the 
temperature  coefficient  of  reactivity  can  be  obtained  experimentally. 
The  remaining  nuclear  coefficients  and  heat  transfer  coefficients  must 
be  calculated  from  the  tabulated  data  for  materials  and  compared  where 
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possible  with  known  physical  characteristics  of  the  reactor  to  es 
lish  their  reliability.  For  examples,  the  nuclear  parameters  should 
predict  the  actual  critical  mass,  and  the  calculated  kinetic  equa- 
tion coefficients,  when  applied  to  a  computer  solution  for  low  power 
transient  behaviour,  should  provide  curves  in  reasonable  agreement  w 
curves  obtained  experimentally. 

The  dependent  variables  in  the  kinetic  equations  _;re  temperature, 
neutron  density,  and  the  concentrations  of  the  six  delayed  neutron 

precursors.  The  effect  of  nuclear  poisons  is  omitted  since  the  peak 

11  2 

flux  at  1,000  watts  is  in  the  order  of  10   neutrons/cm  -sec  at  which 

poison  effects  are  not  significant.   From  neutron  density,  neutron  flux 

and  power  can  be  determined.   In  the  following  sections,  a  description 

of  the  reactor,  the  procedure  for  obtaining  parameters,  and  the  method 

of  obtaining  computer  solutions  for  neutron  density  and  temperature  ai^ 

described. 


2.  Description  of  the  Reactor 

The  reactor  assembly  is  shown  in  Figure  1  and  further  detail  of 
the  core  is  shown  in  Figure  2.   The  core  material  is  a  homogeneous 
mixture  of  polyethylene  and  uranium  dioxide;  the  uranium  is  20  per- 
cent enriched  in  Uranium  235.  The  core  is  surrounded  by  a  high- 
density  graphite  reflector,  a  lead  shield,  and  a  water  shield.   For 
powers  higher  than  100  milliwatts  additional  shielding  is  added  in 
the  form  of  concrete  blocks.   The  structural  material  in  the  core 
region  is  6061-T6  aluminum. 

Control  is  affected  by  the  insertion  and  withdrawal  of  rods  con- 
taining essentially  the  same  material  as  the  core.  Three  rods  are 
the  same  size;  two  of  these  serve  as  safety  rods,  and  the  third  is 
a  coarse  control  rod.   The  fourth  rod,  which  is  smaller,  provides 
fine  control. 

Table  1  summarizes  the  characteristics  of  the  AGN  201  series 
of  reactors  in  general. 
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Figure   1.      AGN  201   Reactor  Assembly 


CORE  TANK  COVER 


GRAPH/ TE  REFLECTOR 


-FUEL  D/SCS 


GLORY  HOLE 


THERMAL  FUSE 


CORE  TANK 


CONTROL  &SAFETY 
ROD  TH/MBLES 

-CORE  SUPPORT ASSEM. 


CORE  TANK  COVER 


Figure  2,   AGN  201  Core  Assembly 


3.  Preliminary  Considerations. 

From  the  characteristics  noted  in  the  previous  section,  it  is  evident 
that  the  reactor  has  little  excess  reactivity  at  room  temperature.   Dur- 
ing laboratory  periods,  it  has  been  found  that  five  or  six  of  the  small 
cadmium  covers  used  around  foils  in  thermal  flux  determinations  are  suf- 
ficient to  shut  the  reactor  down.  This  limited  excess  reactivity  and  the 
strong  negative  temperature  coefficient  would  prevent  continuous  operation 
with  the  glory  hole  empty  at  core  temperatures  above  27  °C.   Since  the  re- 
actor has  no  cooling  system,  the  average  temperature  in  the  core  would  in- 
crease seven  degrees  quite  rapidly  when  the  reactor  is  operating  with  a 
high  neutron  flux  and  the  resulting  high  heat  generation  rate  in  the 
core. 

One  would  expect  then  that  it  is  necessary  to  approach  high  powers 
as  rapidly  as  possible,  that  temperature  would  increase  significantly 
with  a  corresponding  reduction  in  reactivity  when  the  reactor  is  at  suf- 
ficiently elevated  powers,  and  that  in  a  relatively  short  time  a  negative 
reactivity  would  result  causing  a  decrease  in  power.   Subsequent  long 
period  oscillations  of  power  about  an  equilibrium  position  probably  would 
be  the  end  result. 

To  establish  a  set  of  equations  which  predict  this  behaviour  and 
can  be  solved,  it  is  necessary  to  use  a  number  of  approximations.   Fore- 
most of  these  is  the  Diffusion  Equation  approximation  of  the  exact 
Boltzman  Equation;  in  this  approximation  isotropy  of  the  neutron  veloci- 
ties is  assumed.  Also,  exact  energy  distributions  are  replaced  by  the 
appropriate  averages  over  energy  intervals.  Dividing  the  neutron  popu- 
lation into  groups  according  to  the  energy  provides  excellent  results  if 
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a  large  number  of  groups  are  established.  The  complexity  of  solving 
this  problem,  however,  is  not  justified  except  in  the  most  refined  deter- 
minations of  neutron  behaviour  in  materials. 

One-group  theory,  trodified  one-group  theory,  Fermi  Age  theory  and 
two-group  theory  are  frequently  used  as  approximations.  Of  these,  Fermi 
Age  theory  is  a  good  approximation  for  reactors  moderated  by  relatively 
heavy  materials,  while  two-group  theory  provides  reasonable  results  for 
reactors  in  which  hydrogen  is  the  principal  moderating  material.   (1) 
When  the  moderator  contains  hydrogen  and  a  number  of  heavier  materials, 
the  Selengut  -  Goertzel  equation  yields  better  results  than  a  two-group 
approach  since  the  Selengut  -  Goertzel  equation  treats  hydrogen  slowing 
down  and  slowing  down  due  to  other  materials  separately.   However,  to 
obtain  reasonable  parameters  for  the  time  dependent  problem  with  more 
than  one  region,  two-group  theory  is  more  readily  applied. 

A  second  area  of  approximations  is  the  geometric  model  of  the  re- 
actor to  be  assumed.   Since  the  theories  require  the  solution  of  several 
simultaneous  partial  differential  equations  containing  the  Laplacian  and 
two  or  three  independent  variables,  a  solution  for  the  actual  geometry 
would  be  extremely  difficult  to  obtain.   The  problem  is  sufficiently 
complicated  if  one  assumes  a  core  region  surrounded  by  a  reflector  re- 
gion with  cylindrical  geometry,  and  ignores  all  details  of  structural 
materials  and  construction.  This  model  is  satisfactory  from  the  nuclear 
standpoint,  but  further  consideration  is  required  in  connection  with  heat 
transfer,  since  the  structural  aluminum  is  a  good  conductor  and  the 
boundary  and  interface  conditions  have  a  different  physical  significance. 
These  features  are  discussed  in  Section  6. 
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For  a  single  region  containing  one  material,  the  time -dependent 
neutron  diffusion  equation  for  a  single  energy  group  simplifies  to: 

DVH(nt)  -I4>(f,t)  +  SP,t)  =^f±) 

and  for  this  same  region  the  time -dependent  heat  diffusion  equation  is 
given  as : 

KV2T(?,t)  +  ^ftt)-^™ 

In  both  equations,  the  coefficients  may  be  functions  of  the  independent 
variables  in  the  two  equations;  for  example,  the  conductivity  and  dif- 
fusion coefficient  are  both  temperature  dependent.   The  source  term  in 
the  heat  diffusion  equation  is  proportional  to  the  thermal  neutron  flux 
for  the  same  point  in  the  region,  and,  in  turn,  the  source  term  and  the 
coefficients  of  the  neutron  diffusion  equation  are  affected  by  the  tempera- 
ture o 

Considering  a  small  volume  element  in  the  region  using  a  two-group 
approach,  we  have  "the  following  coupled  partial  differential  equations 2 

DfV\-  ^A+  vZMi-/5)  +  %  \Q  =  |?     (1=1  to 6) 


DSV\  -  Lak  +  P  LA  =  ff 


KV2T  +  E<k=Ac 


where  the  subscripts  "f"  and  "s"  refer  to  the  fast  and  slow  neutron 
groups,  respectively,  and  where  E  is  a  factor  to  convert  neutron  flux 
in  the  element  to  the  corresponding  heat  generation  rate   These  equa- 
tions could  be  solved  by  numerical  methods  using  finite  difference 
techniques  and  a  mesh  system.  The  complexity  of  such  a  solution  and 
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Che  time  required  to  solve  this  problem  are  prohibitive  even  with  the  use 
of  a  high  speed  digital  computer.   Further  discussion  of  this  problem  is 
included  in  Section  8, 

A  more  practical  approach,  and  the  approach  used  most  frequently,  is 
to  obtain  from  the  above  equations  a  set  of  ordinary  differential  equa- 
tions amenable  to  computer  solutions.   Considering  only  the  thermal  group 
of  neutrons,  the  diffusion  equation  may  be  written^ 

where  p1  r^  refers  to  the  resonance  escape  probability  and  the  fast 
non-leakage  probability  of  the  delayed  neutrons.  With  the  assumption 
that  the  flux  can  be  calculated  from  the  wave  equation 

V2^(r,t)+B2<t>(M>o 

the  equation  above  may  be  reduced  to 


ft  =  X^"*i£*V:i 


<eff 

(See  Appendix  1) 

This  equation  coupled  with  the  six  ordinary  differential  equations  for 
the  delayed  neutron  precursors  and  a  differential  equation  for  tempera- 
ture can  be  solved  readily  by  numerical  methods. 

Irrespective  of  the  approach,  it  is  necessary  to  obtain  values  for 
coefficients  and  parameters  which  are  required  for  the  equations.   Some 
of  these  parameters  can  be  calculated  quite  accurately,  and  some  data 
are  tabulated.   In  other  cases,  the  values  must  be  estimated.  The  test 
of  any  set  of  parameters  selected  is  that  they  predict  features  of  the 
reactor  which  are  known.   For  example,  the  parameters  should,  when  used 
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In  the  appropriate  equations,  yield  approximately  the  known  critical  mass 
and  the  steady-state  flux  distribution,  and  solutions  of  the  transient 
equations  at  low  power  should  agree  with  corresponding  curves  from  the 
reactor.   For  this  reason  a  number  of  experiments  were  carried  out  on 
the  reactor  to  obtain  data  both  to  supplement  calculated  data  and  to 
provide  verification  of  calculated  data»  These  experiments  are  described 
in  the  following  section,, 
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4,   Experimental  Data  Obtained. 

(a)  Excess  Reactivity  Available 

One  of  the  most  important  parameters  required  for  the  high 
power  estimation  is  the  maximum  excess  reactivity  available,,   Since  this 
quantity  is  very  small,  it  cannot  be  calculated  with  any  degree  of  ac- 
curacy but  may  be  determined  experimentally  with  reasonable  precision  from 
the  reactor  period.   The  rate  at  which  excess  reactivity  can  be  added  is 
also  of  importance.   Two  forms  of  reactivity  can  be  accomplished,  the 
step  input  and  the  ramp  input.   The  ramp  input  is  the  usual  form  of  in- 
put, but  a  step  input  could  be  arranged  if  necessary. 

The  method  of  applying  a  ramp  input  would  be  the  insertion  of  a  rodj, 
and,  since  the  coarse  rod  contributes  more  reactivity  per  unit  length 
than  the  fine  rod,  coarse  rod  movement  only  is  considered.   The  procedure 
then  for  determination  of  excess  reactivity  and  rate  of  application  of  ex- 
cess reactivity  was  as  follows; 

With  the  fine  rod  fully  inserted  and  the  glory  hole  empty,  the  re- 
actor period  was  determined  for  various  positions  of  the  coarse  rod  from 

1 
the  Varian  Recorder  charts.   This  experiment  was  repeated  with  polye- 
thylene in  the  glory  hole.   Finally,  the  actual  rate  of  coarse  rod  travel 
was  established  from  the  average  of  a  number  of  runs  In  which  travel  for 
a  specified  stop  watch  time  was  measured.   The  temperature  of  the  water 

Two  ion  chambers  and  one  BF  counter  are  fitted  in  the  reactor  for 
monitoring  the  power  level.   A  Varian  Recorder  is  used  to  record  the 
micro-micro  ammeter  output  for  either  ion  chamber  or  the  count -rate 
meter  of  the  BF  counter. 
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shield  on  the  day  of  the  experiment  was  (21.2  +  0.5)°C  and  it  is  assumed 
that  the  core  temperature  was  very  nearly  the  same. 

The  rate  of  coarse  rod  movement  was  found  to  be  (26.3  +  0.2)  cm/ 
min. ;  this  figure  coupled  with  an  almost  linear  relation  between  insertion 
and  excess  reactivity  provided  data  for  a  ramp  input.   The  relation  be- 
tween rod  Insertion  and  /~      ,  the  ratio  of  reactivity  to  delayed  neu- 

0 

tron  effectiveness,  is  shown  in  Figure  3(b).   It  is  noted  that  the  contri- 
bution to  reactivity  per  unit  length  of  the  rod  is  the  same  for  polye- 
thylene in  the  glory  hole  as  for  an  empty  glory  hole  for  practical  pur- 
poses.  The  linear  relationship  agrees  with  the  approximately  linear  por- 
tion of  the  standard  S-curve  since  the  rods  of  the  AGN  201  do  not  pass 
through  the  entire  core.   From  the  slope  of  the  curve  in  Figure  3(b), 
the  value  of  ■££-    per  unit  length  is  approximately  5.61  x  10  /cm. 

The  function  of  excess  reactivity   originates  from  the  equation 


A     ~keffT        4-    1  +  XiT 


A. 


Si 


+  !4T^ 


T 


keffT 
The  factor  )$    allows  for  higher  effectiveness  of  the  delayed  neutrons  as 

compared  to  the  prompt  fission  neutrons  in  reaching  thermal  energy.   In 

the  experiments  for  excess  reactivity,  the  periods  used  were  sufficiently 

long  so  that  the  first  term  of  the  equation  was  insignificant.   Since  the 

factor  ft"  is  not  well  knowns  -f~-      was  computed  and  can  be  used  in  this 

0 

form  as  shown  in  Section  7.   Further  discussion  of  the  origin  and  signi- 
ficance of  the  equation  is  given  in  Appendix  1. 

Since  the  worth  per  unit  length  of  the  coarse  rod  is  approximately 
constant  over  the  last  few  centimeters  of  travel,  the  worth  of  the 

14 


polyethylene  in  the  glory  hole  can  be  related  to  a  length  of  coarse  rod. 
Over  a  period  of  time,  the  coarse  rod  position  with  and  without  polye- 
thylene in  the  glory  hole  was  noted  in  connection  with  the  temperature 
coefficient  estimation.   The  difference  of  positions,  the  worth  of  the 
polyethylene  is  3.94  +  ,03  cm, 

(b)  Temperature  Coefficient  of  Reactivity. 

The  temperature  coefficient  experiment  was  actually  a  rough 

-4     1 
attempt  to  check  the  reported  experimental  value  of  (-2.75  x  10  )/°C. 

This  experiment  assumes  that  the  effect  on  the  core  of  a  constant  tempera- 
ture throughout  is  the  same  as  for  the  numerically  equal  average  tempera- 
ture of  a  temperature  distribution  in  the  core.   This  is  not  quite  ac- 
curate because  a  volume  element  near  the  centre  of  the  core  has  a  great- 
er "importance"  than  a  volume  element  near  the  surface  of  the  core  and 
is  at  a  higher  temperature.   If  an  average  temperature  is  used,  it  should 
be  weighted  according  to  the  importance  distribution.   The  experiment  then 
serves  only  to  provide  an  approximation  of  the  magnitude  of  the  tempera- 
ture coefficient  of  reactivity. 

Over  a  period  of  several  weeks,  the  overall  temperature  of  the 
reactor  was  raised  and  lowered  over  as  wide  a  range  as  possible  by  ad- 
justing the  thermostat  of  the  heating  system  in  the  reactor  facility.   The 
temperature  range  achieved  was  only  two  and  a  half  degrees  by  this  elemen- 
tary procedure,  but  sufficient  difference  was  noted  in  the  position  of 
the  coarse  rod  for  criticality  to  make  the  plot  shown  in  Figure  3(a). 
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Figure  3(a)  Temperature  Coefficient  of  Reactivity  Detenn 
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The  least  square  curve  through  the  points  has  a  slope  of  +  .472  cm/°C 

-4 
corresponding  to  a  reactivity  of  -  2.75  x  10  /°C.   This  agrees  sur- 
prisingly well  with  the  value  obtained  by  AGN  for  a  much  wider  range  of 
temperatures. 

At  this  points  the  questions  may  be  raised  as  to  whether  the  worth 
per  centimeter  of  the  coarse  rod  determined  at  21,2°C  is  the  same  for 
other  temperatures  around  20°C  and  if  the  polyethylene  has  an  effect  on 
the  temperature  coefficient  estimate  as  compared  to  the  reactor  with  the 
glory  hole  empty.   A  rerun  of  the  excess  reactivity  experiment  at  a  dif- 
ferent temperature  produced  no  significantly  different  results  and  the 
worth  of  the  polyethylene  over  the  temperature  range  showed  no  consist- 
ent tendency  to  increase  or  decrease  as  a  function  of  temperature. 

(c)  Drop  Test 

A  further  experiment  providing  some  insight  to  the  worth  of 
the  polyethylene  in  the  glory  hole  was  carried  out;  the  primary  purpose 
of  this  experiment  was  to  approximate  a  negative  step  input  so  that  the 
computer  programs  could  be  checked.   This  experiment  consisted  of  obtain- 
ing a  Varian  Recorder  chart  of  the  reactor  behaviour  when  the  polye- 
thylene was  pulled  suddenly  from  the  glory  hole  while  the  reactor  was 
at  a  steady  power.  This  procedure  is  similar  to  the  Drop  Test  method  of 
obtaining  reactivity  for  a  group  of  rods  (1)  (4).   Plotting  the  logar 
of  flux  versus  time  over  the  first  few  seconds  and  extrapolating  back 
from  the  approximate  straight  line  slope  as  shown  in  Figure  4  provides  a 
value  of  flux  q).  ,  on  the  t  -  0  axis.  Then  from  the  equation 
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a  value  for  /— -  of  0.0023  +  0.0002  is  obtained  which  agrees  with  the 

"ft" 

value  0.0022  obtained  from  the  product  of  coarse  rod  worth  per  centimeter 
and  worth  of  polyethylene  in  centimeters  of  coarse  rod. 

(d)   Neutron  Flux  Traverse  through  the  Glory  Hole 

The  relative  thermal  neutron  flux  and  the  epi- cadmium  flux 
through  the  glory  hole  are  shown  in  Figure  5.  The  fluxes  were  obtained 
by  irradiating  bare  and  cadmium- covered  indium  foils  at  various  equi- 
spaced  positions  in  the  glory  hole  and  determining  the  activity  of  the 
radiated  foils  with  a  scintillation  counter.   The  flux  plots  were  consider- 
ed necessary  to  obtain  an  indication  of  geometric  buckling  by  measuring 
the  extrapolated  core  radius.  Also,  from  the  slopes  at  the  interface  and 
the  curvature  in  each  region,  relative  values  of  the  diffusion  coefficients 
and  values  of  diffusion  lengths  could  be  estimated.   The  experimental  pre- 
cision of  these  flux  plots  is  not  particularly  good,  but  the  plots  serve 
the  purposes  for  which  they  were  intended. 

It  is  of  interest  to  note  that  there  is  a  significant  dip  in  the 
thermal  flux  at  the  centre  of  the  core  due  to  the  thermal  fuse  so  that 
the  peak  thermal  flux  is  about  one  centimeter  from  the  center.   Also, 
the  functions  associated  with  the  minor  buckling  of  a  two-group  analysis 
which  are  usually  positive  for  thermal  flux  and  negative  for  fast  flux 
at  the  core-reflector  interface  are  of  reversed  sign  for  the  AGN  201. 
This  fact  was  useful  as  a  guide  to  obtaining  relative  values  of  dif- 
fusion coefficients  and  flux  equation  terms  in  the  parameter  determina- 
tions discussed  in  the  next  section.   One  further  use  made  of  the  flux 
plot  was  a  location  of  the  center  of  the  core;  first  attempts  at  obtain- 
ing a  flux  plot  indicated  that  the  foil  holders  in  use  were  off  by  about 
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one-half  inch  in  calibration  and  subsequent  corrections  were  made. 

The  data  of  this  section  and  the  data  of  the  next  two  sections  are 
summarized  on  page  49  . 
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5.  Theoretical  Nuclear  Parameters. 

As  was  mentioned  previously  in  Section  3,  two-group  theory  pro- 
vides a  reasonable  approximation  for  reactors  moderated  by  hydrogenous 
materials.   If  a  set  of  two-group  parameters  are  obtained  which  predict 
quite  accurately  several  known  features  of  the  reactor,  it  is  assumed 
that  these  parameters  and  other  data  calculated  from  them  are  good  ap- 
proximations for  subsequent  computations. 

Two-group  theory  for  a  reactor  with  a  single  reflector  region  pro- 
vides the  following  diffusion  equations  assuming  no  spatial  dependence 
of  coefficients  within  each  region; 

D,V2<t>,(f)-  Z%(r)  -h  vZ^r)    ■=   0 


Core 


Reflector 


D5V%?)-r(*4vp)  =o 


The  solutions  of  these  equations  provide  equations  for  the  steady- state 
slow  and  fast  flux  at  all  points  in  the  core  and  reflector  and  a  critical- 
ity  determinant  from  which  the  critical  mass  can  be  calculated, 

(a)   Reflector  Parameters 

For  the  reflector,  five  parameters  are  required;   the  thermal 
and  fast  diffusion  coefficients,  the  macroscopic  thermal  cross  section, 
the  resonance  escape  probability  and  the  fast  neutron  removal  cross 
section.   The  removal  cross  section  is  obtained  from  the  ratio  of  the 
fast  diffusion  coefficient  over  the  ficticious  fast  diffusion  length 
(frequently  taken  as  equal  to  the  Fermi  Age). 
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The  reflector  of  the  AGN  201  is  manufactured  from  high-density 

3 
reactor-grade  graphite  with  a  nominal  density  of  1,75  gm/cm  .  A  rough 

3 
check  on  the  density  of  a  sample  yielded  1.74  +  .02  gm/cm  so  that  the 

nominal  figure  is  satisfactory  for  calculations.   The  standard  density 

3 
for  graphite  in  most  references  is  1.60  gm/cm  and  the  tabulated  nuclear 

parameters  for  commerical  graphite  are  based  on  this  density  (1).   Apply- 
ing a  correction  to  the  tabulated  parameters  for  density  differenee9  the 
following  data  are  obtained: 

[)   -  0.838  cm. 

Y*(4)  *  0.000374  cm"1 

* 

£}      m   0.929  cm 

2 
y     s  304  cm 

Then   Y^  =  -E|_  «,  0.00306  cm"1 
The  remaining  factor,  the  resonance  escape  probability   pr s  is  ef- 
fectively 1.00  since  carbon  has  no  resonance  peaks  in  the  epithermal 
range  (5)  and  the  impurity  content  of  reactor-grade  graphite  is  very  low. 

(b)  Core  Parameters  -  First  estimate 

Unfortunately ?  core  parameters  are  not  as  readily  obtained. 
Polyethylene  has  been  used  quite  frequently  in  shielding  and  as  the  modera- 
tor of  critical  assemblies^  but,  for  the  most  part8  the  nuclear  data  for 
water  have  been  used  in  calculations  (with  a  hydrogen  atom  density  ratio 
correction).  Very  little  experimental  nuclear  data  on  polyethylene  it- 
self are  available. 

From  the  official  records  of  "License  Material"  at  the  U.  S.  Naval 

Postgraduate  School,  the  following  data  pertain  to  the  core: 

235 
Total  weight  of  0  665.03  gm 
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235 
Percent  U    In  uranium  19.9370 

235  3 

U    density  in  polyethylene  0.055  gm/cm 

Total  weight  of  components  14781.6  gm 

Total  weight  of  Uranium  3337.67  gms 

The  uranium  is  in  the  form  of  uranium  dioxide  particles  compacted  with 

the  polyethylene. 

235 
The  value  of  665.0  gms  for  the  weight  of  U    is  probably  quite  ac- 
curate, and  for  all  components,  except  the  thermal  fuse,  the  weight  ratio 

235 
of  U    to  total  weight  is  consistent  (0.0450  +  0.0003).   If  the  impregna- 

235 
tion  density  is  0.055  gm  U    cc,  the  overall  density  is  1.22  gm/cc. 

From  the  weight  of  components  and  the  volume  calculated  from  AGN  201 

blueprints,  a  density  of  (1.22  +  0.01)  gm/cc.  is  indicated  for  core 

discs  but  a  slightly  lower  value  is  obtained  for  rod  components.   Since 

the  rods  form  but  a  small  part  of  the  core  and  the  error  involved  Is  less 

than  one  percent,  an  overall  density  of  1.22  gms/cc.   will  be  assumed  for 

the  polyethylene  and  uranium  dioxide. 

3  3 

The  core  also  contains  about  250  cm  of  aluminum,  and  about  50  cm 

of  pure  polyethylene  will  be  in  the  glory  hole  to  obtain  maximum  excess 
reactivity.  The  aluminum  has  little  effect  on  overall  nuclear  character- 
istics of  the  core  since  both  the  absorption  cross  section  and  scattering 
cross  section  of  aluminum  are  relatively  small.   The  polyethylene  contri- 
butes significantly  to  the  excess  reactivity  as  previously  noted,  but  the 
volume  concerned  is  about  0.4%  of  the  core  volume.   The  calculations  of 
the  core  characteristics  can  then  be  based  with  little  error  on  the  uran- 
ium-impregnated polyethylene  with  a  density  of  1.22  gm/cc.   On  this  basis 
the  Table  2  data  are  calculated.   Since  the  volume  ratio  of  U0~  to  the  pol 
ethylene  moderator  is  small,  the  fast  and  thermal  diffusion  coefficients 
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Table  2 
Core  Material  Data 


235 
Weight  of  U 


Core  material  volume 

Overall  density  of  core 

235 
Densities  U 

u238 

Uranium 
Oxygen 

uo2 

Carbon 
Hydrogen 
Polyethylene 
Volume  ratio  of  UP.  to  total  volume  2.97. 


665.03  gm 

12,100  cm3 

1.22  gm/cc 

0.055 

gm/cc 

0.221 

gm/cc 

0.276 

gm/cc 

0.037 

gm/cc 

0.313 

gm/cc 

0.779 

gm/cc 

0.130 

gm/cc 

0.909 

gm/cc 

and  the  Fermi  Age  are  approximately  the  same  as  for  pure  polyethylene. 

Three  sources  of  neutron  cross  sections  are  included  in  the  Biblio- 
graphy and  are  referenced  in  Table  3  below.   Although  the  1960  Supple- 
ment to  BNL  325  is  not  included,  data  from  the  supplement  were  used  in 

reference  (7).   All  microscopic  cross  sections  are  in  bams. 

235 
For  U  *  ,  Westcott's  tables  (6)  will  be  used  since  they  provide  the 

most  up  to  date  data  on  materials  with  non  -  ~—     cross  sections  and  alJow 

for  the  1/E  tail  of  the  neutron  distribution. 

Using  the  densities  of  the  individual  materials  and  the  cross  sections 
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as  listed,  microscopic  absorption  cross  sections  averaged  over  a  Max- 
wellian  distribution  for  a  temperature  of  20°C  are  obtained  from  the 
equation 


E  = 


(«§&)  -tT.L 


1. 128 

Table  3 
Neutron  Cross  sections 


H  (gas 
C 


0 


235 


D 


238 


&« 


c^ 


^ 


<n 


V 


.332(5)    (7) 

38(5)    (7) 

.0034(5) 

4.2(5)    (7) 

.00373(7) 

.0002(5) 

4.2(5)    (7) 

.0002(7) 

694(5) 

10(5)                582(5)          2.07(5) 

2.47(5) 

694(7) 

10(7)                582(7) 

683.04(6) 

577.01(6)   2.0695(6) 

2.4498(6) 

2.71(5)(6)(7) 

.0005(7) 

235 
The  cross  sections  for  U  '   in  the  above  equation  are  obtained  from  West- 

cott's  tables  as 


Westcott's  values  are  those  to  which  tabulated  cross  sections  were 
normalized  and  differ  slightly  from  accepted  values. 
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Of  ■  562.5  bams 

A, 

Oa  s  669.2  barns 
for  r  s  0.03  and  the  S~  cut  -  off.   The  value  of  r  defines  the  pro- 
portion of  epithermal  neutrons  in  the  spectrum;  whereas  this  was  not 
measured  experimentally  for  the  AGN  201,  the  value  is  considered  to 
be  a  reasonable  estimate.   The  S2  cut  -  off  corresponding  to  cutting 
off  the  neutron  spectrum  above  4.95  KT,  is  recommended  in  Westcott's 

tables  for  general  cases. 

=  235  -1 

Then  Za  s  0.0837  cm 

^235"         _i 
2_.p  -  0,0710  cm 

Vs  -1 

and  2_q  ■  0.0231  cm   for  the  oxygen9  carbon 

and  hydrogen  lumped  together. 

With  these  cross  sections  the  product  of  the  thermal  utilization 
and  the  fast  neutrons  produced  by  thermal  fission  per  thermal  absorp- 
tion in  the  fuel  can  be  obtained  as 

^■P  -  1.594 
From  the  equation  (1) 

8 


2T*Y.,    ,  .*l\£ 


r/'1"^  ^ 


i 


Four  figure  accuracy  is  not  justified  but  data  of  the  tables  main- 
tain the  additional  figures. 
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and  values  tabulated  with  the  equation,  an  approximation  can  be  obtain- 
ed of  the  fast  fission  factor  for  the  AGN  201.   As  expected  £   does 
not  differ  significantly  from  1.000.   Finally,  with  the  calculation  of 
the  resonance  escape  probability  "  p  "  using  epithermal  scattering  cross 
sections  and  Figure  14  of  Section  6  -  2  in  Reference  (1),  a  value  of 

"  Kqo    "  can  be  obtained  for  the  core  material. 

pc=  0.960 

kcT^epc35  1-5-30 

The  fast  diffusion  coefficient  for  the  core  material  as  a  function 
of  lethargy  may  be  approximated  using  the  equation 


3|lsl(")(i-A)i 

where  "1"  refers  to  the  nuclide.   Since  Jlc-   varies  with  lethargy,  a 
numerical  integration  was  carried  out  using  BNL  325  curves  and  equal 
lethargy  intervals.   The  contribution  of  oxygen  and  uranium  scattering 
was  neglected  since  the  macroscopic  scattering  cross  sections  of  these 
materials  are  relatively  small.   The  numerical  integration  yielded  a 
value  of  0.623  cm  for  D  . 

To  determine  the  removal  cross  section  2Lr  ,  the  effective  "fast 
diffusion  length"  is  first  obtained.   For  hydrogenous  materials,  the  so- 
called  Fermi  Age  is  approximately  equal  to  the  square  of  the  "fast 
diffusion  length"  (2).   An  experimental  value  for  the  Fermi  Age  of 
polyethylene  has  not  been  published  but  the  Age  of  polyethylene  can  be 
calculated  from  that  of  water  because  the  hydrogen  of  both  accounts  for 
most  of  the  macroscopic  scattering  cross  section  and  the  scattering  cross 
section  of  carbon  and  of  oxygen  are  similar  in  the  epithermal  region. 
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The  principal  difference  between  polyethylene  and  water  is  the  atomic 
density.  The  Fermi  Age  of  a  material  varies  inversely  as  the  product 
of  it:  density  and  its  averaged  microscopic  scattering  cross  section. 
Roughly  then 


?£ 


s 


2 
and  since  "J"  =  27.7  cm  for  water  (2) 


^>  ■  0.920  for  water 

C  =  0.908  for  polyethylene 

—  -l 

2_  -  1.50  cm   for  water 


___  -1 

S   g  -  1.78  cm   for  the  core  material 

2 
then  /T"'  =  19.8  cm  for  polyethylene 

The  removal  cross  section  is  given  by 

2_r  =  -^l  =  0.0237  cm 

For  the  thermal  diffusion  coefficient,  it  is  necessary  to  estimate 
again  from  the  known  data  for  water.   The  thermal  diffusion  coefficient 
for  water  averaged  over  a  thermal  distribution  is  given  as  0.159  cm. 
Correcting  for  the  ratio  of  atomic  densities  provides  a  first  estimate 
of  0.136  cm.  for  polyethylene  if  the  scattering  cross  section  and  the 
average  cosine  of  che  scattering  angles  are  assumed  similar.   However, 
the  scattering  cross  section  of  polyethylene  is  as  much  as  157.  higher 
than  that  of  water  at  thermal  energies  and  exhibits  more  strongly  the 
effects  of  molecular  binding.   If  the  discrepency  was  averaged  over  a 
thermal  distribution  covering  the  entire  thermal  range,  the  difference 
might  be  in  the  order  of  10%  resulting  in  an  estimate  of  0.123  cm  for 
D?.   As  a  check  on  this  supposition,  the  relation 
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3D=  r         ^ 

was  averaged  numerically  over  a  Maxwellian  distribution  with  kT  equal  to 
0.0253  ev.   Linear  relations  were  assumed  for  \\~Mo)    of  water  and  poly- 
ethylene and  2_ t  as  a  function  of  energy  was  taken  from  BNL  325  in  the 
case  of  water  and  from  an  unpublished  KAPL  Paper  for  polyethylene. 
The  ratio  of  the  two  averaged  functions  when  applied  to  the  diffusion 
coefficient  of  water  (0.159  cm)  yielded  a  corresponding  value  for  poly- 
ethylene of  0ol25  cm.   This  value  is  assumed  as  a  reasonable  estimate 
for  purposes  of  a  two-group  analysis. 

(c)  Checking  Parameters  in  Critical  Mass  Calculations. 

With  the  nuclear  parameters  calculated  thus  far  in  this  sec- 
tion, a  two-group  two-region  analysis  for  cylindrical  geometry  was  at- 
tempted for  purposes  of  computing  the  critical  mass  as  one  check  on  the 
validity  of  the  parameters.  The  results  of  the  attempt  were  unsatisfact- 
ory, and  subsequent  attempts  were  made  with  small  variations  in  the  para- 
meters in  an  effort  to  establish  whether  the  poor  results  were  due  to 
minor  inaccuracies  within  the  tolerances  of  the  basic  data. 

The  two-group  analysis  provides  a  four-by-four  determinant,  the 
terms  of  which  consist  of  functions  dependent  on  the  geometry  and  on  the 
coefficients  of  the  partial  differential  equations.   This  determinant  can 
be  expanded  into  an  equation  with  the  most  sensitive  function  of  core 
parameters  on  one  side  of  the  equation  and  the  less  sensitive  functions 
on  the  other  side  (2). 

Private  communication  with  R.  E.  Slovacek,  Knolls  Atomic  Power 
Laboratory„ 
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Making  slight  adjustments  to  various  parameters  of  the  core  in 
either  the  determinant  or  the  equation  to  achieve  a  balance  is  very 
laborious  and  tends  to  obscure  the  significance  of  changes.   The  matrix 
formulation  of  Reference  (1)  however  separates  out  the  parameters  of  the 
reflector,  which  are  fairly  well  known,  into  a  reflector  matrix,  and  the 
adjustment  of  parameters  in  the  core  matrix  is  both  easier  and  more  en- 
lightening. 

When  satisfactory  results  were  still  not  achieved,  the  analysis 
was  repeated  for  spherical  geometry  to  establish  whether  the  "corners" 
omitted  in  the  cylindrical  geometry  were  the  source  of  error.   These 
corners  are  shown  in  Figure  6(a).  The  cases  shown  in  Figure  6(b)  and 
Figure  6(c)  can  be  solved  analytically  and  the  case  in  Figure  6(d)  can 
be  approached  by  iteration  using  the  geometries  of  6(b)  and  6(c)  alter- 
nately.  However,  with  a  small  cylindrical  core  surrounded  by  a  relative- 
ly thick  reflector,  the  corners  omitted  in  6(d)  are  significant.   The  act- 
ual core  dimensions  are  very  close  to  the  optimum  case  for  cylindrical 
geometry,  and  this  in  turn  has  very  similar  characteristics  to  a  sphere. 
If  an  equivalent  spherical  geometry  is  assumed  with  the  same  buckling  as 
the  cylinder  and  approximately  the  same  reflector  savings,  the  omission  of 
a  large  part  of  the  reflector  in  the  cylindrical  geometry  case  can  be 
avoided. 

However s  it  was  found  that  no  justifiable  adjustment  of  the  estimated 
parameters  would  satisfy  the  criticality  condition.   The  parameters  in 
every  case  led  to  a  core  size  considerably  smaller  than  the  actual  core. 

(d)  Effect  of  Spectrum  Hardening 

A  factor  which  had  not  been  included  in  the  calculations  (and 
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Figure  6(a) 


Figure  6(d) 
Figure  6  Geometries  for  Analytical  Solutions 
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ts  not  Included  in  the  calculations  thus  far  for  comparison  purposes) 
was  "neutron  spectrum  hardening".   This  additional  factor,  resolves  the 
problem  of  not  predicting  the  correct  critical  mass. 

Spectrum  hardening  refers  to  a  shift  of  the  curve  for  neutron  dis- 
tribution as  a  function  of  energy  towards  higher  energies.   When  the 
moderating  material  has  a  constant  scattering  cross  section  over  the 
thermal  energy  range,  the  shift  is  due  to  the  -prz-      or  approximately 
absorption  cross  section  of  nuclides  contained.   In  such  cases, 


IT 
the  shift  is  slight  if  the  absorber  concentration  is  not  large.   When 

the  scattering  cross  section  changes  significantly  in  the  thermal  range 

due  to  effects  of  molecular  binding,  the  spectrum  is  further  hardened. 

12  3 
Recent  papers   *  *   indicate  that  for  polyethylene  moderated  assemblies, 

the  shift  is  very  significant  and  leads  to  lower  averaged  macroscopic 

cross  sections.   Neglecting  this  effect  can  result  in  a  large  difference 

in  the  criticality  determinant.   The  spectrum  shift  in  water  can  be  pre- 

4 
dieted  quite  accurately  using  Nelkin's  Kernel  ,  but  no  kernel  has  been 

developed  for  polyethylene. 
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A  survey  of  the  papers  on  this  subject  for  polyethylene-moderat- 
ed assemblies  indicated  that  the  peak  of  the  distribution  should  occur 
at  0.029  or  0.030  electron  volts  for  the  absorber  to  scatterer  ratio 
of  the  AGN  201  core  material  (the  unshifted  spectrum  peak  is  at  0.0253 
electron  volts  for  20°C).   In  calculating  macroscopic  cross  sections 
averaged  over  this  shifted  spectrum,  an  effective  neutron  temperature  of 
336°K  to  347°K  would  be  used.   These  effective  neutron  temperatures 
correspond  to  ratios  of  average  neutron  velocity  to  the  2,200  meter/ sec 
value  of  1.21  to  1.23  as  compared  to  1.128  for  an  unshifted  spectrum. 
The  empirical  formula  (1) 

^  -  1.128  +  1.3G  ^^ 

where  <y     is  the  average  neutron  velocity  of  the  distribution  and  V 
is  the  velocity  corresponding  to  energy  1<0T  yields  a  value  of  1.23 
which  tends  to  confirm  the  data  taken  from  spectrum  plots.   For  purposes 

of  calculations,  assume  for  —  absorbers 

IT 


ct 


1.128  U 


JELtf'  To  -  293°K 

Teff    22°° 


T0 


^OPc\n  where  T   ts   tne  actual 

I  <  ZZ  \J  I    c.ao\j 

temperature 
Applying  this  equation  to  recalculation  of  the  macroscopic  cross  sections 
of  the  core  materials, 

235  .i 


2>a-  .0781 

-^=Z35 

2_f  =  .0646 


cm 


cm 


-  .0214  cm   for  carbon,  oxygen  and 
hydrogen  combined 
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^(f  «  1.590 

Comparing  to  the  previous  calculated  values,  the  effect  on   K  oo  is 
slight,  but  the  absorption  and  fission  macroscopic  cross  sections  are 
decreased  roughly  seven  or  eight  percent.   The  thermal  diffusion  co- 
efficient would  also  be  affected  by  a  spectrum  shift,  but  the  value  for 
water  from  which  it  was  estimated  is  an  experimental  value  and  must  then 
include  this  effect.  The  value  of  0.125  cm  for  the  thermal  diffusion 
coefficient  of  polyethylene  is  at  best  a  guess  in  any  event. 

(e)  Two-Group  Parameters  Which  Predict  Correct  Critical  Mass 

The  two  group  parameters  after  all  corrections  are  summarized 


in  Table  4. 


Table  4 
Two  Group  Nuclear  Parameters 


Dx  m   0.623  cm  Jc   *  19.7  cm       k^  1.52? 

D2  =  0.125  cm  2T^»  0.0995  cm"1    pc  m   0.960 

D  -  0.929  cm  ^f  m  304  cm         pr^  1.00 

D.  -  0.838  cm  T"(4L  0.000374  cm"1 

4  *-a 


2 
Assuming  an  effective  reflector  thickness  of  30  cm  ,  these  parameters 

predict  for  spherical  geometry  a  core  radius  of  14.2  cm.,  a  critical  vol- 

3 
ume  of  approximately  12,000  cm  which  is  just  less  than  the  actual  criti- 
cal volume  and  an  extrapolated  radius  of  approximately  8  inches,  slightly 

For  computational  purposes.,  three  figure  precision  is  assumed,  but  the 
data  does  not  justify  this  degree  of  precision. 

2 
The  graphite  is  20  cm  thick,  but  the  lead  and  water  make  some  contri- 
bution to  neutron  reflection. 
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greater  than  that  indicated  in  the  flux  plot  (which  would  be  expected 
for  an  equivalent  sphere  approximation  ).  The  parameters  then  may  not 
be  quite  correct,  but  are  reasonably  close  to  correct  values. 

(f)  Temperature  Coefficient  of  Reactivity,  Delayed  Neutron  Effect- 
iveness and  Neutron  Lifetime. 

Using  the  parameters  derived  for  the  two-group  analysis, 
several  other  parameters  can  be  estimated.   As  a  check  on  the  experiment- 
al value  of  the  temperature  coefficient  of  reactivity,  an  analytical  value 
may  be  calculated  using  the  method  of  Reference  (2).   From  the  equation 

we  obtain 

J_l!<_       ±VB  +  J_^£  .  !*£.  .  lla.' la.  .  J_llk 

With  the  data  available  and  using  plots  of  various  parameters  as  a 
function  of  temperature  as  required: 

IK  so 

£    3T 
I    3f 


f    3T 


1  m  - L_  MliS1) L_  5f^  _       4  0  x  .^ 

K[    ST  l+Q'     BT  f,+a     ^T  -        4.0*\0  /  c 

t  |f  =  -  -  25.9  x  \0S/oc 
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c     ~    _L  lk_ 

T  k     ^T 

CT    ~    -  2.6  x   I0~4/oc 

Cj      -    -  2.7SX  I0_4/oc 

The  value  of  Cy  compares  c,uite  well  with  the  experimental  value, 
but  depends  heavily  on  the  Fermi  Age  of  the  core  material  which  is  not 
well  known.   Also,  the  theory  is  based  on  an  equivalent  bare  reactor  and 
neglects  any  effects  contributed  by  the  reflector. 

The  overall  effectiveness  of  delayed  neutrons  ^f  can  be  estimated 
from  the  fast  leakage  and  resonance  escape  probability.   Delayed  neutrons 
are  emitted  with  an  average  energy  of  0.515  Mev  (1)  whereas  prompt  neutrons 
have  an  average  energy  of  2  Mev.   The  resonance  range  for  Uranium-238, 
which  is  the  chief  resonance  absorber,  lies  below  0.515  Mev.   The  reson- 
ance escape  probabilities  should  then  be  the  same  for  delayed  and  prompt 
neutrons.   Fast  leakages  for  both  types  of  neutrons  would  also  be  similar 
for  materials  to  which  the  Fermi  Age  theory  may  be  applied,  since  the  age 
is  assumed  to  vary  almost  linearly  with  lethargy.   For  hydrogenous  modera- 
tors the  continuous  slowing  down  model  does  not  apply,  and  the  age  to 
Indium  resonance  in  water  for  an  initial  neutron  energy  of  0.97  Mev„ 
is  roughly  half  the  age  for  the  prompt  neutrons.    Since  0.97  Mev.  is 
representative  of  the  delayed  neutron  initial  energies,  the  ratio  of 
fast  leakage  of  the  delayed  neutrons  to  that  of  the  prompt  neutrons  is 
by  the  two-group  approach  approximately 

D.  Graham  Foster,  Jr.,  Age  of  Na-Be  Neutrons  in  Water  and  Kerosene, 
Nuclear  Science  and  Engineering,  Vol.  8,  No.  2,  Aug.  1960» 
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I  t-  Td  BZ       I  +  1 0.6x0.024 


1/ 


I  +  )fBz       1+  |  0.7x0  024 
The  delayed  neutron  effectiveness  then  is  probably  slightly  greater  than 

1.17  since  the  above  estimate  of  the  age  for  delayed  neutrons  is  too 

La  i  ge. 

One  other  nuclear  parameter  that  will  be  required  is  the  neutron 

lifetime.   For  thermal  neutrons,  this  is  given  approximately  by  the 

formula  , 

A  = 


Slight  corrections  are  necessary  to  include  the  fast  neutron  lifetime  in 
the  overall  neutron  lifetime.   The  method  of  correction  and  the  correc- 
tion used  are  given  in  Section  7. 
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6.   Heat  Transfer  Calculations. 

The  relative  rates  of  heat  generation  in  the  core  and  heat  trans- 
fer from  the  core  at  any  time  set  limits  on  the  future  capability  of  the 
reactor,,   The  heat  generation  rate  is  a  function  of  the  neutron  flux; 
the  heat  transfer  rate  is  a  function  of  temperature.   With  a  high  neu- 
tron flux  and  a  corresponding  high  heat  generation  rate,  the  tempera- 
ture rises  towards  a  value  where  the  heat  transfer  rate  balances  the 
heat  generation  rate.   The  temperature  rise  however  has  an  adverse  effect 
on  nuclear  parameters  such  that  the  ability  to  maintain  or  to  increase 
power  is  reduced  as  the  temperature  increases. 

In  specifying  the  heat  generation  rate  in  the  core,  it  is  assumed 
that  each  fission  of  a  fuel  atom  creates  a  certain  amount  of  energy 
realized  in  the  form  of  heat.   The  value  frequently  used  is  200  Mev.  per 
fission  although  various  values  between  190  Mev.  and  200  Mev.  are  report- 
ed.  About  90%  of  the  energy  is  produced  at  the  time  of  fission;  the  re- 
mainder, resulting  from  the  decay  of  fission  products,  is  produced  at  a 
rate  dependent  on  fission  product  decay  constants.   The  heat  generated 
then  is  not  instantaneous,  nor  is  it  localized.   Depending  on  the  size 
and  construction  of  a  reactor,  a  certain  percentage  of  the  heat  is  pro- 
duced in  each  of  the  fuel,  moderator  and  shield  regions.   For  a  small 
thermal  reactor,  about  90%  of  the  heat  is  produced  in  the  fuel  and  modera- 
tor, the  remainder  in  the  reflector  and  shielding.   (1)   For  computer  soli 
tions,  it  will  be  assumed  that  about  90%  of  the  energy  of  the  fissions 
appears  in  a  very  short  time  in  the  form  of  heat  in  the  core. 

In  an  analysis  of  the  heat  transfer  from  the  core,  effects  such  as 
eddy  currents  in  the  air  gaps  and  the  reduction  of  thickness  of  air  gaps 
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due  to  thermal  expansion  are  negligible.   Over  the  temperature  range 
of  interest,  the  heat  transfer  is  by  conduction  and  the  variation  of 
thermal  conductivity  with  temperature  can  be  disregarded.   These  assump- 
tions lead  to  a  heat  transfer  rate  which  is  a  linear  function  of  the 
temperature  difference  between  the  core  and  whichever  region  remote  from 
the  core  is  assumed  to  remain  at  a  steady  temperature.  . 

The  water  tank  containing  half  a  ton  of  water  could  absorb  all  the 
heat  generated  by  the  reactor  operating  at  1,000  watts  for  one  hour  with 
less  than  a  two  degree  centigrade  temperature  rise.   With  only  short  per- 
iods at  high  power  envisaged,  the  water  tank  serves  as  the  heat  sink  for 
analysis  purposes.   Working  inward  towards  the  core,  the  temperature  dif- 
ference through  the  lead  may  be  estimated  assuming  a  spherical  shell  for 
which 

Using  a  conductivity  for  lead  of  20  BTU/hr-ft-°F,  the  temperature  differ- 
ence is  about  three  degrees  Fahrenheit  per  kilowatt  of  power.   Similarly, 

for  the  graphite,  a  temperature  difference  is  estimated  as  3.7°F  per 

I 
kilowatt  of  power  assuming  a  conductivity  of  109  BTU/sec-ft°F  . 

Next,  since  the  temperature  coefficient  of  reactivity  mentioned  in 

Section  4  is  based  on  the  average  temperature  of  the  core,  the  relation 

between  average  core  temperature  above  the  core  surface  temperature  and 

the  heat  flow  from  the  core  is  required.   A  first  approximation  may  be 


The  conductivity  for  graphite  (with  a  density  of  1.75  gm/cc)  parallel 
to  direction  of  extrusion  is  about  0.55  cal/sec-cm- °C  (1).   Perpendicu- 
lar to  the  direction  of  extrusion,  the  conductivity  is  less  by  30  to 
50%.   The  value  0.45  cal/sec-cm- °C  is  taken  as  a  mean. 
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obtained  assuming  a  spherical  core  with  a  thermal  flux  distribution 
varying  as 
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From  the  equation 


V'TW.^Ifl-^-jf 


r4£r 


br 


sin 


TTT 
Re 


TO)  -  T0-k^rr\z 


\ 


Sin 


Re 


Re 


] 


TTT 
R«* 


Then  the  total  heat  produced  in  the  core  per  unit  time  is 
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and  the  average  temperature  over  the  core  is 
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Relating  the  difference  between  average  and  surface  temperature  to  the 


heat  generation  with  Re  -  20 


cms 


r1  =  14.2  cm  and  using  a  conductivity 


for  polyethylene  of  8x10   cal/sec-cm-°C, 

f-TV, 


a 


0,209  F/Bh.A 


Neglecting  contact  resistances  and  air  gaps  at  this  point,  the  average 

temperature  of  the  core  would  be  over  700°F  greater  than  the  surface 

temperature  if  steady  state  operation  at  1  kilowatt  were  possible.  The 

temperature  drop  across  the  lead  and  graphite  is  therefore  negligible  and 

based  on  this  model,  a  net  thermal  conductance  of  approximately  4.7  BTU/hr- 

°F  may  be  assigned. 
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The  spherical  core  model  does  not  take  Into  account  the  various 
heat  paths  in  the  actual  construction  such  as  the  glory  hole  casing 
and  the  aluminum  control  rod  casings.   Some  indication  of  their  signi- 
ficance can  be  obtained  by  considering  a  network  consisting  of  the  flow 
paths  between  the  various  core  components  and  the  thermal  resistance 
of  the  paths.   Considering  each  disc  of  the  core  separately  and  assum- 
ing a  radial  flux  distribution  proportional  to  the  Bessel  function  J0  » 
the  resistance  to  radial  heat  flow  may  be  obtained  from 

-2-^,       .   t  ,      s  I        ,   2.-404B 


KV  T(r)  -h  ^0J0(<*r)  ro   where  ot» 
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q  =  r  'v  Jo(^od^ 
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Substituting  the  core  radius  r  =  12.9  cm,  the  extrapolated  radius 
R  »  18.8  cm  (from  two-group  analysis)  and  the  conductivity  of  poly- 
ethylene        — 
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From  blueprints  of  the  core  region  of  Figure  2,  there  is  an  air 
gap  between  the  sides  of  the  core  and  the  graphite  approximately  equal 
to  0.1  inches.   At  the  top,  there  is  a  gap  of  approximately  one-quarter 
inch.   The  bottom  disc  rests  on  a  graphite  block.   The  aluminum  thimbles 
and  casings  of  the  control  rods  provide  a  net  area  for  heat  flow  down- 
wards of  about  1.3  square  inches  but  have  air  gaps  around  them  of  ap- 
proximately 0.1  inches  which  tend  to  hinder  heat  flow  to  the  rods. 
Further  heat  paths  are  the  glory  hole  casing  and  the  divider  plate  be- 
tween the  top  and  bottom  half  of  the  core. 

Combining  the  components  and  resistances  between  them,  including 

air  gaps  and  contact  resistance,,  the  network  of  Figure  7  is  obtained. 

2 
Contact  resistances  were  included  on  the  basis  of  h  =  500  BTU/hr-ft  - 

c 

°F.   Each  core  disc  can  lose  heat  to  the  next  disc,  radially  to  the 
graphite,  or  to  an  aluminum  component  if  it  is  near  one,  and  receives 
heat  according  to  its  size  and  to  an  assumed  cosine  axial  distribution 
of  the  heat  generation  rate. 

The  network  can  be  solved  by  the  relaxation  technique  or  a  matrix 
inversion.   The  latter  method  was  chosen,  and  for  an  arbitrarily  chosen 
overall  heat  generation  rate,  the  temperatures  at  each  point  and  the  heat 
flows  along  certain  paths  were  determined.   Results  of  the  analysis  in- 
dicated that  about  707.  of  the  heat  flows  directly  to  the  graphite  out 
the  top,  sides  and  bottom  of  the  core,  the  remainder  through  the  various 
paths  offered  by  the  aluminum  components.  Temperatures  in  the  top  half 
of  the  core  are  higher  than  in  the  bottom  half  due  to  the  lack  of  aluminum 
paths  and  the  poor  conductivity  of  the  air  gap  at  the  top  as  compared  to 
pure  contact  resistance  at  the  bottom.  The  overall  average  temperature 
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Figure  7  Heat  Transfer  Network  for  Core  Components 
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an  assumed  1,000  BUT/hr  total  heat  generation  rate  was  over  100 °F 
corresponding  to  a  conductance  for  the  core  slightly  less  than  10  BTU/hr- 
°F.   This  value  of  conductance  is  about  twice  that  obtained  for  the 
spherical  approximation.   It  is  quite  probable  that  this  approximation 
over-estimates  the  significance  of  the  aluminum  components  and  that  the 
actual  conductance  is  bracketed  by  the  two  approximations. 

Since  it  is  anticipated  that  the  temperature  coefficient  for  a  rapid 
rise  of  power  would  be  greater  than  the  experimentally  determined  value, 
the  lower  estimate  of  conductance  is  favored  for  partial  compensation. 
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7.   Computer  Solution  of  the  Coupled  Differential  Equations, 

The  equation  for  the  neutron  density  at  a  point  in  the  reactor 
if  prompt  and  delayed  neutrons  are  treated  equivalent ly  is: 

(Appendix  1,  Eqn.  1) 
If  the  delayed  neutron  effectiveness  is  included  and  the  assumptions  of 
Appendix  1  are  used,  this  becomes: 

where  /O     is  the  effective  reactivity  at  any  temperature.   With  the 
temperature  dependence  of  the  reactivity  displayed,  this  becomes: 

where  £>       is  the  effective  reactivity  at  an  arbitrary  reference  tempera- 
ture and  AT  is  the  difference  between  the  temperature  of  the  core  and 
the  reference  temperature.  The  corresponding  equations  for  the  six  de- 
layed neutron  group  precursors  have  the  form: 

dCj.  _  _\.r.  +.     ke\+  "X /3i  n 

Jt  ~      c  L  A 

The  temperature  dependence  equation,  assuming  an  average  tempera- 
ture in  the  core  is: 

c1(ct)  =    t*vet  .-  _   KoiAT) 

dt  z3^  Ac 

where  S  is  the  energy  released  per  fission  expressed  in  units  of  heat 
and  K  is  the  conductance  of  the  core  based  on  average  core  temperature 

The  Ra  -  Be  neutron  source  in  the  reactor  is  neglected  since  the 
contribution  to  overall  neutron  density  is  small. 
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(See  Section  6). 

The  eight  coupled  differential  equations  then  have  the  form 

^l-  -X-UC;  4-  G/Si  n 


where 


These  equations  cannot  be  solved  analytically,  and  analog  computer 
solutions  are  limited  to  relatively  small  variations  of  any  of  the 
dependent  variables  unless  special  provisions  are  made.   Using  numerical 
methods,  digital  computer  solutions  for  large  variations  can  be  obtained 
with  any  of  a  number  of  integration  techniques,, 

Of  these  techniques,  the  Euler  method,  modifications  of  Euler6 s 
method  and  the  Runga-Kutta  method  are  frequently  used,,  The  Euler  method 
consists  of  computing  each  new  point  on  the  basis  of  the  last  point  and 
the  slope  of  the  function  evaluated  at  that  point: 


y'-y°-faklsx 


An  improvement  of  Euler8 s  method  consists  of  computing  an  estimated  value 
and  providing  an  improvement  by  averaging  the  slopes  for  last  point  and 
the  estimated  next  point; 


/.  =  Yo  f 


$L)    Ax 


o-> 


"-^iBIU^W 


kx 


0-M 


Both  the  Euler  and  modified  Euler  techniques  lack  the  smoothness  of  the 
Runga-Kutta  technique  which  consists  of  computing  several  estimated  values 
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and  averaging  the  values  according  to  a  numerical  technique  (3).  Other 
methods  provide  greater  improvement,  but  with  more  complexity .   The  Milne 
method  and  the  Hamming  method  base  calculations  for  a  point  on  the  value 
and  slope  of  three  previous  points  as  well  as  applying  a  ''corrector"  and 
"modifier"  formula. 

The  coefficients  of  the  equations  to  be  solved  are  calculated  from 
data  in  the  preceeding  three  sections.   These  data  are  summarized  in 
Table  5.   The  delayed  neutron  parameters  are  listed  in  Table  6.   Using 
a  specific  heat  of  polyethylene  of  0.55  cal/gm-°C  and  a  specific  heat  of 
uranium  dioxide  of  0.057  cal/gm-°C,  the  following  values  are  obtained  for 
the  coefficients  corresponding  to  the  case  of  maximum  reactivity  inputs 


F      =      ^VX       ^        X90%  -  2.31  x  10"7oC/sec-neutron/cm3 

Jsb  =        ^/t(poc)  =   1.0  k  ID"4  sec"1 

^-    =    6-79  *  5-6  X  icf4"  «  3.81  x  io"3 

j^(4)  =  0.43DX  5-6  X  IO"4"  =  2.46  x  10"4/sec 


i\  -  41  x  10   sec 


G   =  ^£X_  .  2.8  x  loSec"1 


Before  attempting  a  solution  with  the  maximum  reactivity  input  values^ 
a  program  was  written  based  on  Euler's  method  for  the  National  Cash 
Register  Company  Model  102-A  digital  computer  to  establish  a  comparison 
between  a  computer  solution  and  the  experimental  results  for  the  small 
negative  step  in  reactivity  mentioned  in  Section  4.   A  conclusion  was 
soon  reached  that  these  equations  could  not  be  solved  satisfactorily  with 
a  relatively  slow  computer.  The  small  time  increments  which  were  used 
required  excessive  computer  time  for  the  calculation  of  points  and  larger 
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Table  5 
Summary  of  Experimental  and  Calculated  Data 


Temperature  Coefficient  of  Reactivity 
Position  of  coarse  rod  at  reference 

temperature  of  20°C  with  polyethylene 
in  glory  hole 
Worth  of  coarse  rod  in  units 

of  reactivity 
Rate  of  coarse  rod  insertion 
Position  of  coarse  rod  when 

fully  inserted 
Delayed  neutron  effectiveness 

factor 
Macroscopic  absorption 

cross  section 
Macroscopic  fission  cross 

section 
Average  slow  neutron  velocity 
Geometric  buckling 
Slow  diffusion  coefficient 
Fast  diffusion  coefficient 

Fermi  Age,  (or  fast  diffusion  length)     X    - 
Conductivity  of  core        4.7  Btu/hr, 
Density  of  polyethylene 
Density  of  Uranium  dioxide 


2.75  x  10  7°C 
17.91  cm  (Fig,  3) 


-4 
5.6  x  10  /cm. 


s  26.3  +  .02  cm/ mm  mi* 

-  24.70  cm. 

-  1.17 


-  0.0781  cm 


-1 


=  0.0646  cm 


-1 


- 

2680  m/sec 

B2 

a 

0.024  cm"1 

D2 

- 

0.125  cm 

Dl 

s 

0.623  cm 

r 

- 

19.7  cm2 

< 

K 

<  10  Btu/hr/ °F 

= 

0.909  gm/cc 

& 

0.313  gm/cc 
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Table  6 

Delayed  Neutron  Parameters 

/3,  s  0.00021  X, =  0.0124  sec" 

/32s  0.00140  \a*  0.0505  sec"1 

/33  -  0.00125  \3^  0.111   sec" 

/34^  0.00253  \4*  0.301   sec"1 

/3S*  0.00074  X5-   1.14  sec"1 

/36^  0.00027  X&s  3.01  sec"1 

time  increments  led  to  instability. 

Since  the  basic  service  library  of  the  Control  Data  Corporation 
1604  digital  computer  contains  a  Runga-Kutta-Gill   subroutine  for  the 
simultaneous  solution  of  coupled  differential  equations^  the  program 

for  a  negative  step  in  reactivity  was  re-written  in  Neliac  Compiler 

2 
language  to  take  advantage  of  the  faster  computer  and  a  smoother  inte- 
gration method. 

A  number  of  runs  to  determine  neutron  density  as  a  function  of  time 

-3 
for  a  step  reactivity  of  -2.21  x  10  '  were  made  with  this  program  to 

3 
determine  suitable  time  increments   for  the  method  and  the  effect  of 

changing  neutron  lifetime  for  which  the  calculated  value  of  41  micro- 
seconds seemed  too  small.  The  results  of  these  runs  are  shown  in  Fig.  8 

The  Runga-Kutta-Gill  method  is  an  adaption  of  the  Runga-Kutta  method 
for  computer  use. 

2 
The  Neliac  Compiler  originated  with  the  Navy  Electronics  Laboratory 

in  San  Diego  and  was  adapted  for  use  on  the  CDC  1604  computer  in  1960, 


3 
Although  increments  for  the  independent  variable  can 

some  cases  to  provide  solutions  within  specified  err 

more  practiced  to  duplicate  a  run  with  1/2  the  selec 

compared  solutions  for  satisfactory  error. 
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Figure  8  Comparison  of  Computer  Solution  with  Experimental 
Results  for  a  Negative  step  in  Reactivity 
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with  a  single  curve,  since  it  was  found  that  time  increments  from  0o001 
sec  to  0,01  sec  provided  the  same  data  to  four  figure  accuracy  and  the 
magnitude  of  neutron  lifetime  could  be  doubled  with  little  effect.  The 
computer  solution  is  almost  exactly  five  per  cent  low  which,  in  view  of 
the  approximations  and  experimental  data  involved,  is  within  acceptable 
limits o  The  discrepancy  is  attributed  to  the  difference  between  an  ex- 
perimentally approximate  step  input  and  the  numerical  method  technique 
application  for  an  actual  step. 

Using  the  same  method  as  for  the  negative  transient  program,  a  new 
program  was  written  to  solve  the  kinetic  equations  for  ramp  and  step  in- 
puts, with  or  without  temperature  dependence,,   Details  of  this  program 
are  given  in  Appendix  2.  The  coefficients  obtained  previously  in  this 
section  when  supplied  to  this  program  lead  to  the  curves  shown  in  Fig, 
9.   (Again  it  was  found  that  the  magnitude  of  neutron  lifetime  had  negli- 
gible effect).  The  neutron  density  corresponding  to  1,000  watts  is  in- 
dicated by  a  horizontal  line  and  it  can  be  seen  that  the  program  indicates 
a  capability  considerably  in  excess  of  this  power.   However,  this  solu- 
tion assumes  that  absorbers,  which  would  reduce  the  available  reactivitys 
are  absent  from  the  glory  hole.  The  power  during  this  transient  period 
remains  above  1,000  watts  for  approximately  two  and  a  half  minutes.   If 
the  reactor  were  held  at  a  thousand  watts  such  that  the  heat  generation 
were  constant,  the  temperature  would  increase  at  a  slower  rate  and  provide 
longer  operating  time. 

The  program  was  re-run  with  a  reactivity  corresponding  to  that  avail- 
able with  the  glory  hole  empty.   The  results  are  shown  in  Fig.  10.   It 
can  be  seen  that  a  longer  time  is  required  to  reach  high  powers  and  that 
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the  temperature  is  higher  when  1,000  watts  is  achieved.  Operating 
time  is  reduced,  but  the  difference  in  reactivity  between  this  case 
and  the  case  where  the  glory  hole  was  filled  with  polyethylene  represents 
a  reasonable  amount  of  absorber  which  can  be  irradiated.   The  reactor 
period  during  the  rise  in  power  agrees  within  a  fraction  of  a  second 
with  the  experimental  value  for  the  same  reactivity. 
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Figure  10  Neutron  Density  and  Temperature  for  Maximum  Re- 
activity ramp  Input. 
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8.   A  Two -Group -Computer  Solution. 

A  computer  solution  to  solve  the  simultaneous  partial  differential 
equations  for  the  reactor  was  written  with  the  intention  of  obtaining 
the  flux  at  every  point  for  checking  parameters  and  to  provide  initial 
conditions  for  a  space  and  time  dependent  transient  solution.   The 
transient  problem, however,  was  never  written  as  it  became  evident  that 
excessive  computer  time  would  be  required  to  achieve  any  results  with 
the  approach  envisaged.   The  reason  for  this  may  be  seen  from  consider- 
ing the  following  finite  difference  equation  example. (3)   From  the  time 
dependent  diffusion  equation 

A  -  ...   —  ~        where  V   is  a  variable 

Sx1   St 

the  finite  difference  equation  is  obtained  as  follows 

If  "S"  exceeds  1/2,  the  presence  of  round  off  or  truncation  errors 
causes  exponentially  increasing  oscillations.  The  value  of  "S"  depends 
on  the  mesh  size  chosen  and  in  the  case  of  cylindrical  geometry  with  a 
radial  and  axial  mesh  spacing  "  A  "  of  one  centimeter  and  a  time  in- 
crement "  A  t"  where 

d*iy.i  M4t\j  ]  -  At  Z^j  *  At  urS+  ^  g 

we  obtain  |    _    A+Z^  —    4Atj£-P     S  O 

La 
the  velocity  in  the  fourth  term  results  in  a  time  interval    t  which 

is  very  small.   Roughly 


• 


4  x  &t  *  220000  XO.I25"     .   . 
I    x  \ 

At  <  3x  icf  6  ^ec 

The  computer  would  have  to  work  through  the  mesh  system  for  each  neu- 
tron group  110,000  times  to  reach  the  first  second  of  simulated  time 
which  would  require  many  minutes  of  computer  time.   Other  approaches 
might  be  practical,  but  to  December  1960,  only  spherical  geometry  pro- 
grams have  been  published. 

The  time -independent  partial  differential  equation  program  was 

1  7 

written  to  use  the  Peaseman  Rackford  method  as  used  in  CURE   .   An 

outline  of  the  program  is  given  in  Appendix  3.   Results  were  not  satis- 
factory in  that  convergence  was  poor  and  instabilities  occured  which 
rapidly  increased  in  magnitude.   For  initial  values,  a  cosine  distribu- 
tion in  the  core  and  a  linear  distribution  in  the  reflector  were  used  and 
the  first  few  passes  tended  to  smooth  this  data  towards  the  experimental 
points  as  shown  in  Figure  11.   However,  during  subsequent  passes,  a  hump 
developed  at  the  core  reflector  interface  and  was  rapidly  amplified.  The 
hump  is  shown  in  Fig.  11,  after  the  program  had  been  operated  ten  times. 


D.  W.  Peaceman  and  H.  H.  Rackford,  Jr.,  "The  Numerical  Solution  of 
Parabolic  and  Elliptic  Differential  Equations  J.   Soc.  Industrial 
and  Applied  Math  3,  1955. 

2 

E.  L.  Wachspress  CURE  A  Generalized  Two  Space  Dimension  Multigroup 

Coding  of  the  704  KAPL  -  1724,  May  1957. 
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9»   Discussion  and  Conclusions 

Much  of  the  data  used  in  this  investigation  lack  precision, 
and  the  methods  used  are  approximations  with  varying  degrees  of  re- 
liability.  In  general,  trends  can  be  predicted  but  absolute  values 
of  neutron  density  and  temperature  are  not  reliable.  Temperature 
values  should  be  within  about  ten  percent  if  the  neutron  density  is 
correctly  predicted,  but  neutron  density  might  be  off  by  a  factor  of 
two  at  high  powers.   Comparison  with  experimental  data,  where  possible, 
was  found  to  be  quite  reasonable,  but  experimental  tests  of  the  high 
power  behaviour  are  prevented  by  the  present  license. 

In  view  of  the  reasonable  agreement  with  experimental  data  where 
checks  were  made,  predictions  using  the  same  basic  data  should  not  be 
too  amiss.   On  this  presumption,  the  computer  solutions  indicate  that, 
using  the  total  reactivity  of  the  rods  with  polyethylene  in  the  glory 
hole,  a  power  in  excess  of  1,000  watts  can  be  obtained  and  that  about 
20,000  watts  is  the  maximum  instantaneous  power  under  any  circumstances. 
Although  a  step  input  could  be  arranged  in  order  to  reach  high  powers  in 
less  time,  a  ramp  input  is  sufficient,  since  the  temperature  does  not 
change  significantly  during  the  period  when  reactivity  is  added.   If  the 
operator  can  settle  the  reactor  at  1,000  watts,  shortly  after  it  achieves 
this  power,  the  total  operating  time  at  this  power  is  estimated  from  the 
equation 


d(&0  __  F      __  KqAT     where   AT  <  — 


Using  the  data  of  Section  7  and  the  temperature  when  1,000  watts  power 
is  achieved  (from  Fig.  9)  as  an  initial  temperature,  the  temperature  at 
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which  the  excess  reactivity  is  overcome  is  achieved  in  approximately 
six  minutes.  The  time  is  not  appreciably  increased  if  a  higher  conduct- 
ivity is  assigned  to  the  core  as  the  rate  of  heat  generation  is  signifi- 
cantly greater  than  the  rate  of  heat  flow  from  the  core.   With  absorbers 
in  the  glory  hole,  the  maximum  reactivity  is  less,  the  temperature  when 
1,000  watts  is  achieved  is  higher,  and  operating  time  is  reduced  as  in- 
dicated by  the  behaviour  shown  in  Figure  10. 

The  calculations  indicate  that  "scramming"  the  reactor  immediately  at 
the  end  of  a  high  power  run  (or  when  power  can  no  longer  be  maintained)  is 
necessary  because  of  the  continued  temperature  increase  that  occurs  if 
the  reactor  is  allowed  to  shut  itself  down  by  temperature  effects.   From 
the  temperature  curve  of  Figure  9,  the  average  temperature  increases  to 
27 °C  above  ambient,  and  the  peak  temperature  approaches  the  softening 
point  of  the  thermal  fuse.   Scramming  the  reactor  reduces  the  heat 
generation  rate  by  a  factor  of  10  in  approximately  13  seconds  and  by  a 
factor  of  100  in  approximately  80  seconds;  this  in  turn  reduces  the  tempera- 
ture buildup  considerably. 

Three  major  fields  for  further  research  are  suggested.   First,  ex- 
tensions and  refinements  of  this  work  would  be  both  informative  and  inter- 
esting, such  as  extending  the  analysis  to  include  spatial  dependence  and 
determining  better  estimates  of  the  temperature  coefficient,  peak  power 
and  flux  and  the  temperature  distribution.   Secondly,  there  are  obvious 
gaps  in  the  knowledge  of  nuclear  parameters  used  in  this  work.   In  parti- 
cular, personnel  of  the  various  laboratories  using  polyethylene  for  criti- 
cal assemblies  need  more  precise  data  on  its  nuclear  characteristics. 
The  thermal  diffusion  coefficient  and  the  Fermi  Age  in  particular  are  not 
established.   Finally,  there  is  a  wide  area  for  research  on  computer 
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solutions  of  the  partial  differential  equations  of  diffusion  and 
heat  transfer  problems. 
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APPENDIX  1 
Reactor  Kinetic  Equations 

The  equation  for  the  neutron  density  is 

where  pl  and  If  are  the  delayed  neutron  resonance  escape  probability 
and  their  leakage  while  slowing  doxm  respectively,,   Assuming  a  solu- 
tion may  be  obtained  from  the  wave  equation 

v2K^t)  +  bN>  1^0  =  0 

Substituting  t  2  _  D/^» 

and  kroEa=  irlfp£ 

^  =  -iaO^B>)<j>  +  (.-/5)k00xapf<)»  i-p'P^x.Ci 

Since  (j)=    HU     ;  ?aq>  (  t-»-L2  Ba)  =r  -~      ; 

When  the  delayed  neutron  effectiveness  is  significantly  different  from 
unity,  considerations  based  on  importance  functions  lead  to  the  equa- 
tion 

at  *  V  J^J        ~X~      *  *XlCt 

(Eqn  1) 
where  K  and  JV  are  functions  of  importance  and  are  for  a  two-group 


1Reactor  Handbook  Vol,,  Physics  AECD  -  3645,  March  1955, 
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analysis  of  the  steady  state: 

k'=i 

Ai         1         +  1 

The  ratio  of    //V   appears  in  the  equation  and  assuming  this  is 
not  too  different  from  the      /j\»  uSuaHy  considered,  the  equation 
is  used  for  reactor  solutions. 

From  the  same  approach  (i.e.  weighting  the  parameter  by  a  factor 
"Q"  )  the  equation 


^+rl 


■eft      ir 


3 1 


I  4-\LT 

is  obtained;  this  equation  was  used  in  the  experimental  reactivity  deter- 
mination.  Since  the  term     A-™         was  negligible, 

7  I  KeW 


4-4 


/Si 


and  Eqn  1  can  be  written  as 

(Eqn  2) 
Since  the  temperature  coefficient  of  reactivity  was  related  experi- 
mentally to  /-r?     ,  the  temperature  dependent  case  may  be  treated  with 


the  equation 

TV 


tlon 


(Eqn  3) 
Using  this  notation,,  the  delayed  neutron  precursor  equations  have  the 
form 

dt  W    ^L 

(Eqn  4) 
The  temperature  depends  on  the  rate  of  heat  generation  /  p  Cj)  £  V 
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and  the  rate  of  heat  leakage.   Therefore: 

c\(/\T)   _    E+vweV  _       Ko  ( AT\ 

clt        ~  /°°°  /°°G  ^      J    (Eq„5) 

the  power  at  any  time  may  be  defined  as 


P)    ^=   __  \  /  (     \ wait- sec  ^ 


All  parameters  except  /O     may  be  treated  as  independent  of  tempera- 
ture with  small  error. 
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APPENDIX  2 
Runga  -  Kutta  Gill  Solution  of  Reactor  Kinetic  Equations 

The  equations  to  be  solved  are  outlined  in  Appendix  1. 
For  a  temperature  dependent  transient,  there  are  eight  simultaneous 
equations  to  be  solved,  for  temperature  independent  cases  there  are  seven,, 

The  program  to  solve  these  equations  was  written  in  Neliac  Compiler 
language  with  two  tapes  required.   The  second  of  these  contains  the  data 
for  delayed  neutron  groups  and  the  operational  program.   The  first  tape 
is  made  up  to  the  users  own  requirements.   Data  which  are  furnished  on 
the  first  tape  include  a  program  starting  address  and  floating  point 
decimal  values  for  the  following  in  order: 

(1)  the  parameter  — -  for  the  reactor  in  question 

(2)  the  maximum  value  of  reactivity  change,  positive  or  negative 

(3)  the  rate  of  change  of  reactivity,  positive  or  negative 
(0  for  step  input) 

(4)  the  initial  average  neutron  density 

K 

(5)  the  conductance  term  -1^2.  (0  for  temperature  independent 

case)  f° 

(6)  temperature  coefficient  of  reactivity 

(7)  heat  generation  constant 

(8)  first  time  increment  for  solution  (  £^"t)) 

(9)  second  time  increment  for  solution  (  /^t^) 

(10)  time  increment  for  print  out  (  /\tp  ) 

(11)  time  at  which  time  increments  for  solutions  are  to  be 
changed  (  "tc  ) 

(12)  time  at  which  run  is  to  end. 
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The  user  must  decide  where  the  program  is  to  be  located  in  memory  and 
what  parameters  to  use  -  then  make  up  a  short  tape  in  compiler  language 
as  indicated.  The  second  tape  remains  the  same.   Two  subroutines  from 
the  basic  library  are  required,  the  Runga-Kutta-Gill  subroutine  (loca- 
ted at  60,200)  and  "Decof",  the  floating  point  output  subroutine  (located 
at  71,000).  The  location  of  these  routines  and  the  compiler  location 
place  limits  on  where  the  programs  may  be  located,  but  the  compiler  can 
be  relocated  if  necessary. 

In  operation,  the  program  dumps  on  magnetic  tape  the  values  fed  in 
on  the  first  tape  for  reference  purposes,  then  computes  initial  values 

for  temperature  and  the  delayed  neutron  percursors,  establishes  the  condi- 

K 
tions  for  a  Runga  Kutta  Gill  solution  and,  based  on  whether  1—2.  is  zero 

or  has  a  value,  selects  the  routine  for  a  temperature  independent  or  a 
temperature  dependent  problem.   Until  time  "tc  ,  it  uses  the  time  in- 
crement &-t,  ,  then  changes  to  use  the  time  increment  ^to  •  This  feature 
was  included  such  that  a  shorter  time  increment  for  the  first  few  seconds 
of  a  step  input  problem  could  be  used.   After  each  time  increment  2.*/^tp 
four  or  six  values  are  dumped  on  magnetic  tape.   For  the  temperature  in- 
dependent case,  two  consecutive  times  (in  seconds)  and  the  correspond- 
ing neutron  densities  are  dumped,  for  the  temperature  dependent  case, 
the  temperature  at  each  time  is  also  dumped. 
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APPENDIX  3 
A  Two-Group,  Two-Region  Computer  Program  with  Cylindrical  Geometry 

For  each  region  and  group,  the  neutron  behaviour  in  steady  state 
conditions  is  defined  by  equation  of  the  form 

DV"  K(r,z)  -2Z<f>L(r5z)  +■  50(r>z)  =  0 

where  the  coefficients  are  assumed  constant  throughout  the  region,  and 
the  equations  in  a  region  are  coupled  through  the  source  terms   ^0(r5z) 
where  "o"  signifies  "other  group". 

These  equations  may  be  replaced  by  finite  difference  equations  of 
the  form  (3) 

where  the  coefficients  are  functions  of  position,  mesh  spacing  and  the 
coefficients  of  the  partial  differential  equations. 

For  the  reactor  problem,  a  core  diameter  and  height  of  12.5  cm 

was  selected  with  an  extrapolated  diameter  and  height  for  the  shielding 

JK  z 
of  42.5  cm.   These  dimensions  where 

chosen  to  allow  one  centimetre  spac-  \ 

I 
ing  with  zero  mesh  points  on  the  out-  | 

side  boundaries  and  on  the  core- 
reflector  interface.   Due  to 
symmetry,  only  1/4  of  the  con- 
figuration is  required.   The 
total  mesh  consists  of  44  rows 
and  44  columns  with  1936  mesh  points.  aSvp 

i   44  and  column  44  along  the  ex 
values,  row  1  and  column  1  beside  inte 
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image  values  across  the  axis. 

The  method  of  solution  used  was  the  alternating-direction  technique 

1 
outlined  by  Peaceman  &  Rackford,    In  this  method,  all  the  rows  are 

selected  in  turn  and  solved  by  the  Gaussian  elimination  technique  for 
linear  equations  with  a  coefficient  matrix  of  triple  diagonal  form. 
The  procedure  is  repeated  for  the  columns  in  turn.   Solving  the  rows 
and  columns  once  constitutes  one  pass  of  the  routine.   During  the  solu- 
tion for  any  one  row  or  column,  the  values  in  the  adjacent  row  or  column 
respectively  are  assumed  known  such  that  the  five  point  difference  equa- 
tion has  the  form 

where  "m"  signifies  the  row  (or  column)  which  is  to  be  solved,  don\y)  an<* 
tfU^^are  functions  of  the  coefficients  Qnvyl,  bn   >    ,  Y\    and  dy-j  m 
and  an  iteration  parameter  which  is  varied  depending  on  which  pass  is 
being  executed  (the  parameter  is  large  for  initial  passes,  lower  as  the 
number  of  passes  increases). 

The  program  for  this  problem  was  written  in  Neliac  Compiler  language 
and  consisted  of 

(1)  Storage  areas  for  coefficients,  computed  data  and  intei- 
mediate  data. 

(2)  Routine  for  calculating  all  of  the  coefficients  a   ,  b  , 

nm   nm 

etc.,  and  storing  for  future  use. 

(3)  A  Peaceman- Rackford  method  subroutine,  entered  from  the  main 


in  and  H.  H.  Rachford,  Jr.,  "The  Numerical  Solution 
I  Elliptic  Differential  Equations",  J,  Soc.  Indust. 
3(1955),  p.  28. 
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(4)  An  "Equation"  subroutine  to  solve  a  row  or  column, 
entered  from  the  Peaceman-Rackford  method  subroutine. 

(5)  Main  program,  which  alternately  selected  computation 

of  fast  group  and  slow  group,  provided  checks  and  dumped 
computer  data. 

Such  a  program  would  not  be  expected  to  be  short  due  to  the  storage 
and  extensive  transfer  requirements.  Approximately  12,000  cells  were 
required  for  storage  and  a  further  2,200  foi'  routines  of  which  most 
were  required  for  the  coefficient  calculation  and  the  Feaceman  - 
Rackford  method  subroutine.   One  more  versed  with  computer  programming 
could  reduce  these  requirements  considerably.   What  was  lost  in  the  space 
sense  was  gained  in  time  by  having  many  coefficients  in  storage  "on  call" 
rather  than  calculating  them  as  required.  A  complete  pass  through  the 
Peaceman  -  Rackford  method  routine,  which  involves  calculation  and  trans- 
fer on  two  blocks  of  1,939  values  each  using  the  equation  subroutine, 
required  about  4  seconds. 

In  operation,  convergence  difficulties  were  encountered  which  could 
not  be  eliminated  by  extensive  "de  -  bugging"  and  may  indicate  the  nec- 
essity of  a  new  approach. 
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